Algorytm LU faktoryzacji
Algorytm LU faktoryzacji został wynaleziony przez polskiego matematyka Tadeusza Banachiewicza w roku 1938, co opisano w pracy [1]
Zauważmy iż algorytm eliminacji Gaussa generuje macierz trójkątną górną \( U \). Możliwe jest również wygenerowanie macierzy trójkątnej dolnej \( L \) takiej że oryginalna macierz naszego problemu \( A=LU \). Jaka jest zaleta w posiadaniu takiej dekompozycji? Załóżmy że chcemy rozwiązać dwa problemy dla dwóch prawych stron \( b_1, b_2 \). Załóżmy również że na początku znamy prawą stronę \( b_1 \) a że drugą prawą stronę \( b_2 \) poznamy dopiero po rozwiązaniu problemu dla pierwszej prawej strony. Czy da się w jakiś sposób wykorzystać faktoryzację macierzy \( A \) wykonaną dla poprzedniej prawej strony \( b_1 \)? Najpierw obliczamy LU faktoryzację macierzy \( A=LU \) następnię mamy równanie \( LUx=b_1 \). Obliczamy najpierw \( z \) wykonując podstawianie w przód na podstawie równania \( Lz=b_1 \). Następnie obliczamy x rozwiązując równanie \( Ux=z \) przez podstawianie wstecz.
Zobaczmy jak zmodyfikować algorytm eliminacji Gaussa tak żeby generował obie macierze \( L \) i \( U \). W celu uzyskania macierzy \( U \) uruchamiamy zwykły algorytm eliminacji Gaussa (ale bez prawej strony). Macierz \( L \) uzyskujemy poprzez zapamiętanie liczb przez które mnożymy poszczególne wiersze podczas eliminacji.
Początkowo
\( L= \begin{bmatrix} 1 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 \\ \end{bmatrix} \)
oraz
\( U= \begin{bmatrix} 0.05 & 0.10833 & 0.00833 & 0 & 0 \\ 0.10833 & 0.5 & 0.21666 & 0.00833 & 0 \\ 0.00833 & 0.21666 & 0.55 & 0.21666 & 0.00833 \\ 0 & 0.00833 & 0.21666 & 0.5 & 0.10833 \\ 0 & 0 & 0.00833 & 0.10833 & 0.05\\ \end{bmatrix} \)
Podczas każdego mnożenia aktualizujemy macierz \( L \)
\( 1^{st} = 1^{st} \frac{1}{A_{1,1}}= 1^{st} \frac{1}{0.05}=1^{st}20 \)
\( U = \begin{bmatrix} {\color{blue}0.05 }*20 & 0.10833*20 & 0.00833*20 & 0*20 & 0*20 \\ 0.10833 & 0.5 & 0.21666 & 0.00833 & 0 \\ 0.00833 & 0.21666 & 0.55 & 0.21666 & 0.00833 \\ 0 & 0.00833 & 0.21666 & 0.5 & 0.10833 \\ 0 & 0 & 0.00833 & 0.10833 & 0.05 \\ \end{bmatrix} \)
W szczególności na przekątną wstawiamy wyrazy z przekątnej oryginalnej macierzy.
\( L= \begin{bmatrix} {\color{blue}0.05} & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 \\ \end{bmatrix} \)
\( U= \begin{bmatrix} 1 & 2.16666 & 0.16666 & 0 & 0 \\ {\color{blue}0.10833} & 0.5 & 0.21666 & 0.00833 & 0 \\ 0.00833 & 0.21666 & 0.55 & 0.21666 & 0.00833 \\ 0 & 0.00833 & 0.21666 & 0.5 & 0.10833 \\ 0 & 0 & 0.00833 & 0.10833 & 0.05\\ \end{bmatrix} \)
Podczas odejmowania wierszy, do macierzy \( U \) wstawiamy liczbę używaną do skalowania wiersza
\( 2^{nd} = 2^{nd} - 1^{st} *A_{2,1}= 2^{nd} - 1^{st} *{\color{blue}0.10833} \).
\( L= \begin{bmatrix} 0.05 & 0 & 0 & 0 & 0 \\ 0.10833 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 \\ \end{bmatrix} \)
\( U \begin{bmatrix} 1 & 2.16666 & 0.16666 & 0 & 0 \\ 0.10833-1*0.10833 & 0.5-2.16666*0.10833 & 0.21666-0.16666*0.10833 & 0.00833 & 0 \\ 0.00833 & 0.21666 & 0.55 & 0.21666 & 0.00833 \\ 0 & 0.00833 & 0.21666 & 0.5 & 0.10833 \\ 0 & 0 & 0.00833 & 0.10833 & 0.05\\ \end{bmatrix} \)
\( U= \begin{bmatrix} 1 & 2.16666 & 0.16666 & 0 & 0 \\ 0 & 0.26528 & 0.19860 & 0.00833 & 0 \\ {\color{blue}0.00833} & 0.21666 & 0.55 & 0.21666 & 0.00833 \\ 0 & 0.00833 & 0.21666 & 0.5 & 0.10833 \\ 0 & 0 & 0.00833 & 0.10833 & 0.05 \\ \end{bmatrix} \)
\( 3^{rd} = 3^{rd} - 1^{st} A_{3,1}= 3^{rd} - 1^{st} {\color{blue}0.00833} \)
\( L= \begin{bmatrix} 0.05 & 0 & 0 & 0 & 0 \\ 0.10833 & 1 & 0 & 0 & 0 \\ {\color{blue}0.00833} & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 \\ \end{bmatrix} \)
\( U= \begin{bmatrix} 1 & 2.16666 & 0.16666 & 0 & 0 \\ 0 & 0.26528 & 0.19860 & 0.00833 & 0 \\ 0.00833-1*0.00833 & 0.21666-2.16666*0.00833 & 0.55-0.16666*0.00833 & 0.21666 & 0.00833 \\ 0 & 0.00833 & 0.21666 & 0.5 & 0.10833 \\ 0 & 0 & 0.00833 & 0.10833 & 0.05 \\ \end{bmatrix} \)
\( U= \begin{bmatrix} 1 & 2.16666 & 0.16666 & 0 & 0 \\ 0 & 0.26528 & 0.19860 & 0.00833 & 0 \\ & 0.19861 & 0.54861 & 0.21666 & 0.00833 \\ 0 & 0.00833 & 0.21666 & 0.5 & 0.10833\\ 0 & 0 & 0.00833 & 0.10833 & 0.05\\ \end{bmatrix} \)
\( 2^{nd} = 2^{nd} \frac{1}{A_{2,2}}=2^{nd}\frac{1}{0.26528}=2^{nd}3.76963 \)
\( L= \begin{bmatrix} 0.05 & 0 & 0 & 0 & 0 \\ 0.10833 & {\color{blue}0.26528} & 0 & 0 & 0 \\ 0.00833 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 \\ \end{bmatrix} \)
\( U= \begin{bmatrix} 1 & 2.16666 & 0.16666 & 0 & 0 \\ 0 & {\color{blue}0.26528}*3.76963 & 0.19860*3.76963 & 0.00833*3.76963 & 0 \\ 0 & 0.19861 & 0.54861 & 0.21666 & 0.00833 \\ 0 & 0.00833 & 0.21666 & 0.5 & 0.10833 \\ 0 & 0 & 0.00833 & 0.10833 & 0.05 \\ \end{bmatrix} \)
\( U= \begin{bmatrix} 1 & 2.16666 & 0.16666 & 0 & 0 \\ 0 & 1 & 0.74864 & 0.03140 & 0 \\ 0 & {\color{blue}0.19861} & 0.54861 & 0.21666 & 0.00833 \\ 0 & 0.00833 & 0.21666 & 0.5 & 0.10833\\ 0 & 0 & 0.00833 & 0.10833 & 0.05 \\ \end{bmatrix} \)
\( 3^{rd} = 3^{rd} - 2^{nd} A_{3,2}= 3^{rd} - 2^{st}{\color{blue}0.19861} \)
\( L= \begin{bmatrix} 0.05 & 0 & 0 & 0 & 0 \\ 0.10833 & 0.26528 & 0 & 0 & 0 \\ 0.00833 & {\color{blue}0.19861} & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 \\ \end{bmatrix} \)
\( U= \begin{bmatrix} 1 & 2.16666 & 0.16666 & 0 & 0 \\ 0 & 1 & 0.74864 & 0.03140 & 0 \\ 0 & 0.19861-1*0.19861 & 0.54861-0.74864*0.19861 & 0.21666-0.03140*0.19861 & 0.00833 \\ 0 & 0.00833 & 0.21666 & 0.5 & 0.10833 \\ 0 & 0 & 0.00833 & 0.10833 & 0.05 \\ \end{bmatrix} \)
\( U= \begin{bmatrix} 1 & 2.16666 & 0.16666 & 0 & 0 \\ 0 & 1 & 0.74864 & 0.03140 & 0 \\ 0 & 0 & 0.39992 & 0.21042 & 0.00833 \\ 0 & {\color{blue}0.00833} & 0.21666 & 0.5 & 0.10833 \\ 0 & 0 & 0.00833 & 0.10833 & 0.05 \\ \end{bmatrix} \)
\( 4^{th}=4^{th}-2^{nd}A_{4,2}=4^{th}-2^{nd}{\color{blue}0.00833} \)
\( L= \begin{bmatrix} 0.05 & 0 & 0 & 0 & 0 \\ 0.10833 & 0.26528 & 0 & 0 & 0 \\ 0.00833 & 0.19861 & 1 & 0 & 0 \\ 0 & {\color{blue}0.00833} & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 \\ \end{bmatrix} \)
\( U= \begin{bmatrix} 1 & 2.16666 & 0.16666 & 0 & 0 \\ 0 & 1 & 0.74864 & 0.03140 & 0 \\ 0 & 0 & 0.39992 & 0.21042 & 0.00833 \\ 0 & 0.00833-1*0.00833 & 0.21666-0.74864*0.00833 & 0.5-0.03140*0.00833 & 0.10833 \\ 0 & 0 & 0.00833 & 0.10833 & 0.05 \\ \end{bmatrix} \)
\( U= \begin{bmatrix} 1 & 2.16666 & 0.16666 & 0 & 0 \\ 0 & 1 & 0.74864 & 0.03140 & 0 \\ 0 & 0 & 0.39992 & 0.21042 & 0.00833 \\ 0 & 0 & 0.21042 & 0.499738 &0.10833 \\ 0 & 0 & 0.00833 & 0.10833 & 0.05 \\ \end{bmatrix} \)
\( 3^{rd}=3^{rd}\frac{1}{A_{3,3}}=3^{rd}\frac{1}{0.39992}=3^{rd}2.50050 \)
\( L=\begin{bmatrix} 0.05 & 0 & 0 & 0 & 0 \\ 0.10833 & 0.26528 & 0 & 0 & 0 \\ 0.00833 & 0.19861 & {\color{blue}0.39992} & 0 & 0 \\ 0 & 0.00833 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 \\ \end{bmatrix} \)
\( U= \begin{bmatrix} 1 & 2.16666 & 0.16666 & 0 & 0 \\ 0 & 1 & 0.74864 & 0.03140 & 0 \\ 0 & 0 & {\color{blue}0.39992}*2.50050 & 0.21042 & 0.00833*2.50050 \\ 0 & 0 & 0.21042 & 0.499738 & 0.10833 \\ 0 & 0 & 0.00833 & 0.10833 & 0.05 \\ \end{bmatrix} \)
\( U= \begin{bmatrix} 1 & 2.16666 & 0.16666 & 0 & 0 \\ 0 & 1 & 0.74864 & 0.03140 & 0 \\ 0 & 0 & 1 & 0.52615 & 0.02082 \\ 0 & 0 & {\color{blue}0.21042} & 0.499738 & 0.10833 \\ 0 & 0 & 0.00833 & 0.10833 & 0.05 \\ \end{bmatrix} \)
\( 4^{rd}=4^{rd}-3^{rd}A_{4,3}=4^{rd}-3^{rd}{\color{blue}0.21042} \)
\( L= \begin{bmatrix} 0.05 & 0 & 0 & 0 & 0 \\ 0.10833 & 0.26528 & 0 & 0 & 0 \\ 0.00833 & 0.19861 & 0.39992 & 0 & 0 \\ 0 & 0.00833 & {\color{blue}0.21042} & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 \\ \end{bmatrix} \)
\( U= \begin{bmatrix} 1 & 2.16666 & 0.16666 & 0 & 0 \\ 0 & 1 & 0.74864 & 0.03140 & 0 \\ 0 & 0 & 1 & 0.52615 & 0.02082 \\ 0 & 0 & 0.21042-1*0.21042 & 0.499738-0.52615*0.21042 & 0.10833-0.02082*0.21042 \\ 0 & 0 & 0.00833 & 0.10833 & 0.05 \\ \end{bmatrix} \)
\( U= \begin{bmatrix} 1 & 2.16666 & 0.16666 & 0 & 0 \\ 0 & 1 & 0.74864 & 0.03140 & 0 \\ 0 & 0 & 1 & 0.52615 & 0.02082 \\ 0 & 0 & 0 & 0.38902 & 0.10395 \\ 0 & 0 & {\color{blue}0.00833} & 0.10833 & 0.05 \\ \end{bmatrix} \)
\( 5^{th}=5^{th}-3^{rd}*A_{5,3}=5^{th}-3^{rd}{\color{blue}0.00833} \)
\( L= \begin{bmatrix} 0.05 & 0 & 0 & 0 & 0 \\ 0.10833 & 0.26528 & 0 & 0 & 0 \\ 0.00833 & 0.19861 & 0.39992 & 0 & 0 \\ 0 & 0.00833 & 0.21042 & 1 & 0 \\ 0 & 0 & {\color{blue}0.00833} & 0 & 1 \\ \end{bmatrix} \)
\( U= \begin{bmatrix} 1 & 2.16666 & 0.16666 & 0 & 0 \\ 0 & 1 & 0.74864 & 0.03140 & 0 \\ 0 & 0 & 1 & 0.52615 & 0.02082 \\ 0 & 0 & 0 & 0.38902 & 0.10395 \\ 0 & 0 & 0.00833-1*0.0.00833 & 0.10833-0.52615*0.00833 & 0.05-0.02082*0.00833 \\ \end{bmatrix} \)
\( U= \begin{bmatrix} 1 & 2.16666 & 0.16666 & 0 & 0 \\ 0 & 1 & 0.74864 & 0.03140 & 0 \\ 0 & 0 & 1 & 0.52615 & 0.02082 \\ 0 & 0 & 0 & 0.38902 & 0.10395 \\ 0 & 0 & 0 & 0.10394 & 0.04982 \\ \end{bmatrix} \)
\( 4^{th}=4^{th}\frac{1}{A_{4,4}}=4^{th}\frac{1}{0.38902}=4^{th}2.57056 \)
\( L= \begin{bmatrix} 0.05 & 0 & 0 & 0 & 0 \\ 0.10833 & 0.26528 & 0 & 0 & 0 \\ 0.00833 & 0.19861 & 0.39992 & 0 & 0 \\ 0 & 0.00833 & 0.21042 & {\color{blue}0.38902} & 0 \\ 0 & 0 & 0.00833 & 0 & 1 \\ \end{bmatrix} \)
\( U= \begin{bmatrix} 1 & 2.16666 & 0.16666 & 0 & 0 \\ 0 & 1 & 0.74864 & 0.03140 & 0 \\ 0 & 0 & 1 & 0.52615 & 0.02082 \\ 0 & 0 & 0 & {\color{blue}2.57056}*0.38902 & 2.57056*0.10395 \\ 0 & 0 & 0 & 0.10394 & 0,04982 \\ \end{bmatrix} \)
\( U= \begin{bmatrix} 1 & 2.16666 & 0.16666 & 0 & 0 \\ 0 & 1 & 0.74864 & 0.03140 & 0 \\ 0 & 0 & 1 & 0.52615 & 0.02082 \\ 0 & 0 & 0 & 1 & 0.26721 \\ 0 & 0 & 0 & {\color{blue}0.10394} & 0.04982 \\ \end{bmatrix} \)
\( 5^{th}=5^{th}-4^{th}*{\color{blue}0.10394} \)
\( L= \begin{bmatrix} 0.05 & 0 & 0 & 0 & 0 \\ 0.10833 & 0.26528 & 0 & 0 & 0 \\ 0.00833 & 0.19861 & 0.39992 & 0 & 0 \\ 0 & 0.00833 & 0.21042 & 0.38902 & 0 \\ 0 & 0 & 0.00833 & {\color{blue}0.10394} & 1 \\ \end{bmatrix} \)
\( U= \begin{bmatrix} 1 & 2.16666 & 0.16666 & 0 & 0 \\ 0 & 1 & 0.74864 & 0.03140 & 0 \\ 0 & 0 & 1 & 0.52615 & 0.02082 \\ 0 & 0 & 0 & 1 & 0.26721 \\ 0 & 0 & 0 & 0.010394-1*0.010394 & 0.04982-0.26721*0.010394 \\ \end{bmatrix} \)
\( U= \begin{bmatrix} 1 & 2.16666 & 0.16666 & 0 & 0 \\ 0 & 1 & 0.74864 & 0.03140 & 0 \\ 0 & 0 & 1 & 0.52615 & 0.02082 \\ 0 & 0 & 0 & 1 & 0.26721 \\ 0 & 0 & 0 & 0 & 0.04704 \\ \end{bmatrix} \)
\( 5^{th}=5^{th}-\frac{1}{A_{5,5}}=5^{th}\frac{1}{0,04704}=5^{th}*21.2585 \)
\( L= \begin{bmatrix} 0.05 & 0 & 0 & 0 & 0 \\ 0.10833 & 0.26528 & 0 & 0 & 0 \\ 0.00833 & 0.19861 & 0.39992 & 0 & 0 \\ 0 & 0.00833 & 0.21042 & 0.38902 & 0 \\ 0 & 0 & 0.00833 & 0.10394 & {\color{blue}0.047042} \\ \end{bmatrix} \)
\( U= \begin{bmatrix} 1 & 2.16666 & 0.16666 & 0 & 0 \\ 0 & 1 & 0.74864 & 0.03140 & 0 \\ 0 & 0 & 1 & 0.52615 & 0.02082 \\ 0 & 0 & 0 & 1 & 0.26721 \\ 0 & 0 & 0 & 0\ & {\color{blue}0.047042}*21.2585 \\ \end{bmatrix} \)
W efekcie dostaliśmy następujące macierze \( L \) i \( U \).
\( L * U = \begin{bmatrix} 0.05 & 0 & 0 & 0 & 0 \\ 0.10833 & 0.26528 & 0 & 0 & 0 \\ 0.00833 & 0.19861 & 0.39992 & 0 & 0 \\ 0 & 0.00833 & 0.21042 & 0.38902 & 0 \\ 0 & 0 & 0.00833 & 0.10394 & 0.047042 \\ \end{bmatrix} * \begin{bmatrix} 1 & 2.16666 & 0.16666 & 0 & 0 \\ 0 & 1 & 0.74864 & 0.03140 & 0 \\ 0 & 0 & 1 & 0.52615 & 0.02082 \\ 0 & 0 & 0 & 1 & 0.26721 \\ 0 & 0 & 0 & 0\ & 1 \\ \end{bmatrix} \)
Żeby zweryfikować nasz algorytm, możemy sprawdzić na koniec że \( LU=A \).